///////////////////////////////////////////////////////////////////////////////
//   DKL2ROOT  -  to convert Dackel binary files to ROOT-files
//
//   base version: September 2006 Alex Kaplan
//   ROOT-version: September 2006 Marc Schumann
//   for information on the ZIP-Library see http://www.zlib.net/
//
//   last modifications:
//     13/09/06 MS - main program written
//     14/09/06 MS - Valid flag added to tree
//		   - Pedestal subtraction
//     20/09/06 MS - new Dackel header style included
//		   - modus/StatusCode is read from header and stored in tree
//		   - Counter Mode implemented: append counter values to file
//     03/03/07 TS - change to make easy time cuts possible
//     25/02/07 AK - dackel ev_count bug workarround
//     28/02/07 AK - dkl2root has problems with files containing many events
//     22/03/07 AK - several changes on the number of TDC_STOPS
//     11/04/07 BM - integrate seperate dkl2spec (ASCII spectra) program
//                 - remove useless parameter "r" (CreateASCII CreateROOT)
//                 - merge CreateAscii(), CreateRoot() CreateSpec() into Convert()
//                 - dkl2root now compilable without ROOT (new define HAVE_ROOT)
//     12/04/07 BM - skip invalid cycles, display warning (old and new dackel files)
//                 - new option -z to set the ROOT compression level
//                 - added Root Tree for cycle numbers
//                 - allow substraction of pedestals for spectrum output
//                 - added Root branches for number of trigger events (PmtTrigger),
//                   latch trigger info (Det1Triggered/Det1Triggered,
//                   Det1TriggeredFirst/Det2TriggeredFirst
//
//    -> further history in SVN log ...
//
//   TODO:
//        - add TDC hit counter Tree (attention: different versions!)
//        - handle files with fileOK = 1
//        - add option to reverse TDC info
//        - add options to select SpinUp/Spindown, Det1First/Det2First, Det1/Det2  (ASCII/Spec)
//        - add options to create DeltaTimerCounter cuts
//        - change Valid tree data type to boolean, add some meaning ;-)
//
///////////////////////////////////////////////////////////////////////////////
// Format of ROOT-Tree:
//  |--------------------------------------------------------------------------
//  |Name of Leaf	Size	  Content
//  |--------------------------------------------------------------------------
//  |StatusCode		short	  Status of Cycle coded bitwise (see below)
//  |Valid		short	  Is cycle valid for analysis? Always 1 at the moment!
//  |Latch		integer	  Latch content
//  |TimerCounter	integer	  Time stamp of event. ("occurrence in cycle")
//  |Det1[0..5]		short	  The 6 PMT signals of detector 1
//  |Det2[0..5]		short	  The 6 PMT signals of detector 2
//  |Det1Sum/Det2Sum float   The sum of all 6 PMTs of each detector (scaled)
//  |Det1Tdc/Det2Tdc integer The TDC time of each detector (for backscattering)
//  |--------------------------------------------------------------------------
//
//  Attention: The mapping of the ADCs to the Detector PMTs is done in
//             Function "GetAdcContent()"
//
////////////////////////////////////////////////////////////////////////////////


#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#ifdef WIN32
#  include <io.h>
#  include "getopt.h"
#else
#  include <unistd.h>
#  include <getopt.h>
#endif
#include <string.h>
#include <fcntl.h>
#include <zlib.h>
#include <iostream>
#include <fstream>
#include <cmath>

#ifdef HAVE_ROOT
/* ROOT libraries */
# include <TFile.h>
# include <TH1.h>
# include <TLine.h>
# include <TF1.h>
# include <TROOT.h>
# include <TCanvas.h>
# include <TStyle.h>
# include <TApplication.h>
# include <TGaxis.h>
# include <TTree.h>
# include <TStopwatch.h>
#endif

#include "dackel_date.h"
#include "DklHeader.h"
#include "EventStruct.h"
#include "EventStructVer1.h"
#include "counters.h"

#define DET_NUMBER	2
#define PMT_PER_DET	6

//added 20070322 by alex
#define TDC_STOPS 4
#include "tdc_channels.h"

// TODO: for testing purposes only!!
//#define DET1_CHANNEL 0
//#define DET2_CHANNEL 1

//added 20070402 by alex
// (valid_data & INVALID_DATA) should be zero
// (invalid_data & INVALID_DATA) differs from zero
// see manual for CAEN V767A user's manual!
#define INVALID_DATA 0x00600000

#define PEDESTAL_MAX	1024
#define PEDESTAL_REBIN	1
#define PEDESTAL_RANGE 	20
// minimum time difference between last event and pedestal event
#define PEDESTAL_DELTATIME 30

using namespace std;


#ifndef HAVE_ROOT
// we define some types we use ourselves, since we don't have Root
typedef unsigned int UInt_t;
typedef short int Short_t;
typedef unsigned short int UShort_t;
typedef bool Bool_t;
typedef float Float_t;
#endif


// globale Variablen
char *in_fn = 0;	// input filename
char *out_fn = 0;	// ASCII output filename
char *spec_fn = 0;  // ROOT output filename
char *root_fn = 0;  // ROOT output filename
char *ped_fn = 0;	// Pedestal log-filename
char *counter_fn = 0;	// counter output-filename

#ifdef HAVE_ROOT
// ROOT compression level
unsigned int compressionLevel =    1;  // 1: apply compression to ROOT file
#endif

int in;
gzFile gz_in;
#ifdef HAVE_ROOT
int SubtractPedestals = 0;
#endif
Short_t Pedestal[DET_NUMBER][PMT_PER_DET];
unsigned int ev_count = 0;

// Funktionsprototypen
void show_help(char *progname);
void proc_opts(int argc,char *argv[]);

u_int16_t GetAdcContent(EventStruct *es, int det, int pm, int version);
u_int32_t GetTcContent(EventStruct *es, int version);
u_int32_t GetLatchContent(EventStruct *es, int version);
void GetTdcContent(EventStruct *es, u_int32_t det, UInt_t *dest, int version);
void ReadEvent(EventStruct *es, int version, int i);

#ifdef HAVE_ROOT
int DeterminePedestals(DklHeader &header);
#endif
int Convert(DklHeader &header);
//int CreateAscii(DklHeader header);
//int CreateRoot(DklHeader header);
//int CreateSpec(DklHeader header, char* out_fn);
int WriteCounters(DklHeader &header);


#ifndef HAVE_ROUND
inline double round(double d)
{
  return floor(d + 0.5);
}
#endif


/*****************************************************************************/
int main(int argc,char *argv[])
/*****************************************************************************/
{

  printf("----------------------------------------------------\n"
        "dkl2root by AK, MS, BM    built %s %s\n\n",  __DATE__, __TIME__);

/*
  printf("DEBUG MSG: sizeof(EventStruct)=%u sizeof(EventStructVer1)=%u\n",
        sizeof(EventStruct),sizeof(EventStructVer1));
*/

  // Analysiere Programmparameter
  proc_opts(argc,argv);

  //are there any arguments left?
  if(optind==argc) {
    printf("ERROR: incorrect number of parameters, try %s --help\n\n",
          argv[0]);
    exit(1);
  }

  // get input file name
  in_fn = (char*)malloc(sizeof(char)*(strlen(argv[optind])+1));
  if(!in_fn) {
    printf("error: malloc returned NULL (file:%s line:%d)\n\n",
          __FILE__, __LINE__ );
    exit(1);
  }
  strcpy(in_fn,argv[optind]);

#ifdef HAVE_ROOT
  // check if pedestal is to be substracted and create log-file name
  if (SubtractPedestals && !ped_fn) {
    ped_fn = (char *) malloc(sizeof(char) * (strlen(in_fn) + 5));
    if (!ped_fn) {
      printf("error: malloc returned NULL (file:%s line:%d)\n\n",
             __FILE__, __LINE__ );
      exit(1);
    }
    strcpy(ped_fn, in_fn);
    strcat(ped_fn, ".ped");
  }
#endif

  //open input file
  if( (in=open(in_fn,O_RDONLY | O_BINARY | O_SEQUENTIAL)) == -1 ) {
      printf("error: could not open input file \"%s\"\n\n",in_fn);
      exit(1);
  }

  //read dkl_header
  DklHeader header(in);

  //print info from header
  //    cerr << header.infoString();

  if(header.ev_count == 0 ) {
      cout << "error: event counter == 0\n";
      exit(1);
  }

  // Get Counter Values and put them on screen
  cout << "Cycle\t";
  for (u_int32_t j=0; j<header.noc; j++) cout << "Cnt" << j << "\t"; cout << endl;
  for (u_int32_t i=0; i<header.cycle_count; i++) {
    cout << header.fcID+i << "\t";
    for (u_int32_t j=0; j<header.noc; j++)
      cout << header.getCounter(j,i) << "\t";
    cout << endl;
  }
  if (counter_fn) WriteCounters(header);

  unsigned int r;
  if(	(r=lseek(in,header.hdr_length,SEEK_SET)) != header.hdr_length ) {
    printf("error: lseek returned %d instead of %d\n\n", r, header.hdr_length);
  }

  // initialisiere Pedestals mit Null
  for (int det=0; det<DET_NUMBER; det++)
    for (int pm=0; pm<PMT_PER_DET; pm++)
      Pedestal[det][pm] = 0;
#ifdef HAVE_ROOT
  // bestimme Pedestals
  if (ped_fn) DeterminePedestals(header);
#endif

  if ((out_fn!=NULL) || (spec_fn!=NULL) || (root_fn!=NULL))
    Convert(header);

  //close input file
  close(in);

  printf("\n\n");
  return 0;
}
/*****************************************************************************/


/*---------------------------------------------------------------------------*/
void show_help(char *progname)
/*---------------------------------------------------------------------------*/
{
  printf("Program to convert Dackel binaries to ASCII or ROOT-files\n"
         "usage: %s [options] <inputfile>\n\noptions:\n"
         " -h, --help        shows this help text\n"
         " -o, --data-file   create ASCII file (*)\n"
#ifdef HAVE_ROOT
         " -r, --root-file   create ROOT file (*)\n"
#endif
         " -s, --spec-file   create ASCII spectrum file (*)\n"
#ifdef HAVE_ROOT
         " -p, --pedestals   determine and subtract pedestals\n"
         "                   (a logfile is created automatically)\n"
         " -P                like -p but save data to given file name (*)\n"
#endif
         " -c, --counter     save counter values (*)\n"
#ifdef HAVE_ROOT
         " -z,  --compression-level\n"
         "                   set ROOT compression level\n"
#endif
         "\n(*) options requires a filename as argument\n", progname);
}
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
// Analysiere Uebergabeparameter
void proc_opts(int largc, char *largv[])
/*---------------------------------------------------------------------------*/
{
#ifdef HAVE_GETOPT_LONG
  static struct option long_options[] =
  {
     {"help",              no_argument,       0, 'h'},
     {"data-file",         required_argument, 0, 'o'},
     {"root-file",         required_argument, 0, 'r'},
     {"spec-file",         required_argument, 0, 's'},
     {"pedestals",         no_argument,       0, 'p'},
     {"compression-level", required_argument, 0, 'z'},
     {"counter",           required_argument, 0, 'c'},
     {0, 0, 0, 0}  // muss am Ende immer stehen!
  };
#endif

  optind=0; 	// important! if getopt_long() is used to processes
  		// more than one array, as done in main()
  while(1) {
    /* getopt_long stores the option index here. */
    int option_index = 0;
    int c;

    // important: without this statement getopt_long() does not
    //            terminate processing of arguments correctly
    if(optind>=largc) break;

#ifdef HAVE_GETOPT_LONG
#  ifdef HAVE_ROOT
    c = getopt_long(largc, largv, "ho:r:s:c:pz:P:", long_options, &option_index);
#  else
    c = getopt_long(largc, largv, "ho:s:c:", long_options, &option_index);
#  endif
#else
#  ifdef HAVE_ROOT
    c = getopt(largc, largv, "ho:r:s:c:pz:P:");
#  else
    c = getopt(largc, largv, "ho:s:c:");
#  endif
#endif

    if (c == -1) break; 	// Detect the end of the options.
				// this does not work correctly if you use
				// cr.fargc and cr.fargv! (see above)

    switch(c)
    {
      case 'h': // help-file
         show_help(largv[0]);
	       exit(0);
	       break;
      case 'o': // ASCII-Output
	       out_fn = (char *)malloc(strlen(optarg)+1);
	       if(out_fn) strcpy(out_fn,optarg);
	       break;
#ifdef HAVE_ROOT
      case 'r': // ROOT-Output
        root_fn = (char *)malloc(strlen(optarg)+1);
 	      if(root_fn) strcpy(root_fn,optarg);
	       break;
#endif
      case 's': // ASCII spectrum
        spec_fn = (char *)malloc(strlen(optarg)+1);
 	      if (spec_fn) strcpy(spec_fn, optarg);
	      break;
#ifdef HAVE_ROOT
      case 'P':
        ped_fn = strdup(optarg);
      case 'p': // calculate Pedestals
        SubtractPedestals = 1;
	      break;
#endif
      case 'c': // save counter values
        counter_fn = (char *)malloc(strlen(optarg)+1);
 	      if(counter_fn) strcpy(counter_fn,optarg);
	      break;
#ifdef HAVE_ROOT
      case 'z': // set ROOT compression level
        compressionLevel = atoi(optarg);
        break;
#endif
      case '?':
      	exit(1);
	      break;
      default:
	break;
    }
  } // end while(1)
}
/*---------------------------------------------------------------------------*/


// Gebe ADC-Wert aus Event zurück
/*---------------------------------------------------------------------------*/
u_int16_t GetAdcContent(EventStruct *es, int det, int pm, int version)
/*---------------------------------------------------------------------------*/
{
  // TODO: this function is a mess - please clean up ;-)
  /*
#if 1
  // Bastians C-style version of the function. Kind of a hack, though ;-)
  u_int16_t *p;
  p = &(es.adc00) + 6*det + pm;
  return *p;
#else
  if (det==0) {
  if (pm==0) return es.adc00;
  if (pm==1) return es.adc01;
  if (pm==2) return es.adc02;
  if (pm==3) return es.adc10;
  if (pm==4) return es.adc11;
  if (pm==5) return es.adc12;
  }
  else if (det==1) {
    if (pm==0) return es.adc20;
    if (pm==1) return es.adc21;
    if (pm==2) return es.adc22;
    if (pm==3) return es.adc30;
    if (pm==4) return es.adc31;
    if (pm==5) return es.adc32;
  }
  cout << "error in GetAdcContent(): Det" << det << "Pm" << pm << "does not exist!" << endl;
  exit(1);
#endif
 */
  if(version == 2 ) {
    if (det==0 && pm==0) return es->adc00;
    if (det==0 && pm==1) return es->adc01;
    if (det==0 && pm==2) return es->adc02;
    if (det==0 && pm==3) return es->adc10;
    if (det==0 && pm==4) return es->adc11;
    if (det==0 && pm==5) return es->adc12;
    if (det==1 && pm==0) return es->adc20;
    if (det==1 && pm==1) return es->adc21;
    if (det==1 && pm==2) return es->adc22;
    if (det==1 && pm==3) return es->adc30;
    if (det==1 && pm==4) return es->adc31;
    if (det==1 && pm==5) return es->adc32;

    cout << "error in GetAdcContent(): Det" << det << "Pm" << pm << "does not exist!" << endl;
    exit(1);
  }
  if(version == 1 ) {
    EventStructVer1 *es1 = (EventStructVer1*)es;

    if (det==0 && pm==0) return es1->adc00;
    if (det==0 && pm==1) return es1->adc01;
    if (det==0 && pm==2) return es1->adc02;
    if (det==0 && pm==3) return es1->adc10;
    if (det==0 && pm==4) return es1->adc11;
    if (det==0 && pm==5) return es1->adc12;
    if (det==1 && pm==0) return es1->adc20;
    if (det==1 && pm==1) return es1->adc21;
    if (det==1 && pm==2) return es1->adc22;
    if (det==1 && pm==3) return es1->adc30;
    if (det==1 && pm==4) return es1->adc31;
    if (det==1 && pm==5) return es1->adc32;

    cout << "error in GetAdcContent(): Det" << det << "Pm" << pm << "does not exist!" << endl;
    exit(1);
  }

  printf("GetAdcContent(): unsupported dkl version!");
  exit(1);
  return -1;
}
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
u_int32_t GetTcContent(EventStruct *es, int version)
/*---------------------------------------------------------------------------*/
{
		if (version == 2)
				return es->tc;

		if (version == 1) {
				EventStructVer1 *es1 = (EventStructVer1*) es;
				return es1->tc;
		}

		printf("GetTcContent(): unsupported dkl version!");
		return -1;
}
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
u_int32_t GetLatchContent(EventStruct *es, int version)
/*---------------------------------------------------------------------------*/
{
		if (version == 2)
				return es->latch;

		if (version == 1) {
				EventStructVer1 *es1 = (EventStructVer1*) es;
				return es1->latch;
		}

		printf("unsupported dkl version!");
		return -1;
}
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
void GetTdcContent(EventStruct *es, u_int32_t channel, UInt_t *dest, int version)
/*---------------------------------------------------------------------------*/
{
  if(version == 2 ) {
    u_int32_t *tdc = &(es->tdc1);
    int j=0;

    for(int i=0; i<TDC_STOPS; i++) {
      u_int32_t ch = ((tdc[i] & 0x7f000000) >> 24);
      if( (INVALID_DATA & tdc[i]) ) {
         cerr << "GetTdcContent(): error: invalid datum:"
              << hex << tdc[i] << endl;
      }

      if( ch == channel) {
        dest[j] = tdc[i] & 0xfffff;
        j++;
      }
    }
    //fill rest with zeros
    for (;j<TDC_STOPS;j++)
      dest[j] = 0;

  } else if (version == 1) {
    EventStructVer1 *es1 = (EventStructVer1*)es;
    u_int32_t *tdc = &(es1->tdc1);

    int j=0;
    for(int i=0; i<2; i++) { // in version 1 only 2 TDC_STOPS were saved
      u_int32_t ch = ((tdc[i] & 0x7f000000) >> 24);
      if ((INVALID_DATA & tdc[i])) {
        cerr << "GetTdcContent(): error: invalid datum:"
             << hex << tdc[i] << endl;
      } else if (ch == channel) {
        dest[j] = tdc[i] & 0xfffff;
        j++;
      }
    }
    //fill rest with zeros
    for(;j<TDC_STOPS;j++)
      dest[j] = 0;

  } else {
    printf("GetTdcContent(): unsupported dkl version!");
    exit(1);
  }
}
/*---------------------------------------------------------------------------*/


#ifdef HAVE_ROOT
// Bestimme die Pedestals der ADCs
/*---------------------------------------------------------------------------*/
int DeterminePedestals(DklHeader &header)
/*---------------------------------------------------------------------------*/
/*
  Global variables:
     ped_fn   :   filename for the outout of the pedestal values (logfile)
*/
{
  UInt_t   Latch;
  UInt_t   Cycle;
  unsigned int cyclecounter = 0;
  unsigned int consistency = 0;
  unsigned int lastValidCycle = 0;
  u_int32_t prevLatch = 0;
  UInt_t   TimerCounter;
  UInt_t   DeltaTimerCounter = 0;

  cout << "Calculate Pedestals" << endl;

  // create histograms for pedestal-fits
  TH1D *hist[DET_NUMBER][PMT_PER_DET];
  char HistName[15];

  for (int det=0; det < DET_NUMBER; det++) {
    for (int pm=0; pm < PMT_PER_DET; pm++) {
      sprintf(HistName, "HistDet%uPm%u", det, pm);
      hist[det][pm] = new TH1D(HistName, HistName, PEDESTAL_MAX / PEDESTAL_REBIN, 0, PEDESTAL_MAX);
      if (hist[det][pm] == NULL) {
        cerr << "ERROR: Calculate Pedestals. Cannot allocate memory for histogram Det" << det
             << " PM" << pm << endl;
        exit(1);
      }
      hist[det][pm]->Sumw2();
    }
  }

  // oeffne log-file
  FILE *out;
  if (!(out = fopen(ped_fn, "w"))){
    printf("error: could not open pedestal-log-file \"%s\"\n\n", ped_fn);
    exit(1);
  }

  if (!(gz_in = gzdopen(dup(in), "rb"))) {
    printf("error: gzopen() returned NULL\n\n");
    exit(1);
  }

  EventStruct *es = new EventStruct;
  u_int32_t i=0,p=header.ev_count/49;

#ifdef HAVE_ANSI
  printf("read data for pedestals:\n  [%51s\n\033[1A[","]");		// Fortschrittsbalken
#else
  printf("read data for pedestals:\n  [%51s","]");		// Fortschrittsbalken
  for (i=0; i<51; i++) printf("\b");
#endif

  Cycle = header.fcID - 1;
  cyclecounter--;

  // FIXME: this is not meaningful?!
  lastValidCycle = header.cycle_count-1;
  if (header.fileOK == 0) {
    lastValidCycle = header.cycle_count-1;
  }

  // Fuelle Histogramme
  //for(i=0; i < header.ev_count; i++) {
  i = 0;
  while(!gzeof(gz_in)) {

    if ((i % p) == 0) {
      printf("*");
      fflush(stdout);
    }
    i++;

    ReadEvent(es, header.version, i);
    Latch = GetLatchContent(es, header.version);

	TimerCounter = GetTcContent(es, header.version);
    if (DeltaTimerCounter < TimerCounter)
		DeltaTimerCounter = TimerCounter - DeltaTimerCounter;
    else
        DeltaTimerCounter = TimerCounter;

    // detect new cycles
    if ((Latch >> 31) != (prevLatch >> 31)) {
      cyclecounter++;  // index starting from zero
      Cycle++;         // absolute number (offset header.fcID)
      consistency = 0;
      if (header.noc > CYCLE_VALID_COUNTER) {
        // only newer file versions support this!
        // consistency should be zero for valid cycles, non-zero otherwise
        consistency = header.getCounter(CYCLE_VALID_COUNTER, cyclecounter);
      }
      // skip invalid cycles from aborted files (old dackel versions)
      if (cyclecounter > lastValidCycle)
        consistency = 1;  // TODO: set a more meaningfull value here!
      // report skipped cycles
      if (consistency)
        cout << "Warning: cycle " << Cycle << " (index " << cyclecounter << ") has errors, skipping" << endl;
    }

	if ((!consistency) && (DeltaTimerCounter > PEDESTAL_DELTATIME)) {
		for (int det=0; det < DET_NUMBER; det++) {
			UInt_t   DetTdc[TDC_STOPS];
			if (det == 0)
				GetTdcContent(es, DET1_CHANNEL, DetTdc, header.version);
			else
				GetTdcContent(es, DET2_CHANNEL, DetTdc, header.version);

			if (DetTdc[0]==0)
				for (int pm=0; pm < PMT_PER_DET; pm++)
					hist[det][pm]->Fill(GetAdcContent(es, det, pm, header.version));
		}
	}

    // remember last latch value
    prevLatch = Latch;
  }

  // Fuehre Gauss-Fits an Pedestal durch
  TF1 *gauss;
  char outstring[100];

  for (int det = 0; det < DET_NUMBER; det++) {
    for (int pm = 0; pm < PMT_PER_DET; pm++) {
      int max = (int) hist[det][pm]->GetMaximumBin() * PEDESTAL_REBIN;
      if (max != PEDESTAL_REBIN) {
        hist[det][pm]->Fit("gaus", "", "", max - PEDESTAL_RANGE, max + PEDESTAL_RANGE);
        gauss = hist[det][pm]->GetFunction("gaus");
        sprintf(outstring,"Det%dPm%d: Ped = %.2f +- %.2f    rChi^2 = %.2f/%d = %.2f\n",
                det,pm,gauss->GetParameter(1), gauss->GetParError(1),
                gauss->GetChisquare(), gauss->GetNDF(), gauss->GetChisquare() / gauss->GetNDF());
        Pedestal[det][pm] = (Short_t) round(gauss->GetParameter(1));
      } else {
        sprintf(outstring,"Det%dPm%d: WARNING: no pedestal found; set pedestal to 0\n",det,pm);
        Pedestal[det][pm] = 0;
      }
      cout << outstring;
      fprintf(out, "%s", outstring);
    }
  }

  for (int det=0; det<DET_NUMBER; det++) {
    for (int pm = 0; pm<PMT_PER_DET; pm++) {
      cout << "Pedestal[" << det<< "," << pm <<"] = " << Pedestal[det][pm] << endl;
    }
  }
  fclose(out);

  return 0;
}
/*---------------------------------------------------------------------------*/
#endif


/*---------------------------------------------------------------------------*/
void ReadEvent(EventStruct *es, int version, int i)
/*---------------------------------------------------------------------------*/
{
  int r;

  if (version == 2) {
      if((r = gzread(gz_in, es, sizeof(EventStruct))) <= 0) {
        printf("error: gzread() returned %d, read %u events\n"
                "  gztell() returns %u\n\n", r, i,
                (unsigned int) gztell(gz_in));
        exit(1);
      }
      if(r != sizeof(EventStruct)) {
        printf("warning: gzread() returned %u instead of %u\n"
                " gztell() returns %u\n\n", r, sizeof(es),
                (unsigned int) gztell(gz_in) );
      }
  } else if (version == 1) {
      if ((r = gzread(gz_in, es, sizeof(EventStructVer1)))<=0) {
        printf("error: gzread() returned %d, read %u events\n"
                "  gztell() returns %u\n\n", r,i,
                (unsigned int) gztell(gz_in));
        exit(1);
      }
      if (r != sizeof(EventStructVer1)) {
        printf("warning: gzread() returned %u instead of %u\n"
                " gztell() returns %u\n\n", r, sizeof(es),
                (unsigned int) gztell(gz_in) );
      }
  } else {
      printf("ReadEvent(): Error: dkl file version is %ui, but max supported dkl file version is %ui!", version, 2);
      exit(1);
  }
}
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
int Convert(DklHeader &header)
/*---------------------------------------------------------------------------*/
{
  EventStruct *es = new EventStruct;
  u_int32_t i = 0;
  u_int32_t p = header.ev_count/49;
  unsigned int cyclecounter = 0;
  unsigned int consistency = 0;
  unsigned int lastValidCycle = 0;
  u_int32_t prevLatch = 0;

  UShort_t StatusCode;
  Bool_t   Valid;
  UInt_t   TimerCounter;
  UInt_t   Latch;
  UInt_t   DeltaTimerCounter = 0;
  UInt_t   Cycle;
  UInt_t   TdcHits;
  //Short_t OldDet10 = 0;
  Bool_t   SpinFlipState;
  Short_t  Pmt[DET_NUMBER][PMT_PER_DET];  // must be signed to allow pedestal subtraction !
  Float_t  DetSum[DET_NUMBER];
  UInt_t   Det1Tdc[TDC_STOPS];
  UInt_t   Det2Tdc[TDC_STOPS];
  Bool_t   DetTriggered[DET_NUMBER];
  Bool_t   DetTriggeredFirst[DET_NUMBER];
  Short_t  PmtTrigger[DET_NUMBER];


  // ASCII output
  FILE *out = NULL;

#ifdef HAVE_ROOT
  // variables needed only for ROOT output
  // Variablen, um den Tree zu fuellen
  TFile *out_root = NULL;
  TTree *tree;
  TStopwatch timer;
  // unsigned int maxTreeSize = 100000000;  // maximum size in bytes of the ROOT file
  // unsigned int autoSave    =  50000000;  // write changes automatically every autoSave bytes
#endif

  // SPEC  variables
  FILE *out_spec = NULL;
  //TODO: make this adjustable parameters
  const unsigned spec_size = 65536;
  const unsigned spec_ofs  = 1024;
  const unsigned binsize = 1;  // 1 means no rebinning ;-)
  // allocate spectra
  unsigned pm[2][6][spec_size+1];
  unsigned sum1[spec_size+1], sum2[spec_size+1], sum[spec_size+1];


#ifdef HAVE_ROOT
  // start timer clock
  timer.Start(1);
#endif

  // oeffne Dkl-Datei, wenn nicht schon fuer Pedestal geschehen
  if (!(ped_fn)) {
    if (!(gz_in = gzdopen(dup(in), "rb"))) {
      printf("error: gzopen() returned NULL\n\n");
      exit(1);
    }
  } else { // ... sonst springe nur zurück an Anfang
    if (gzrewind(gz_in) == -1) {
      printf("error in gzrewind()\n\n");
      exit(1);
    }
  }

  // create ASCII output file
  if (out_fn) {
    if (!(out = fopen(out_fn, "w"))) {
      printf("error: could not open output file \"%s\"\n\n", out_fn);
      exit(1);
    }
    fprintf(out, "# pm11 pm12 pm13 pm14 pm15 pm16 pm21 pm22 pm23 pm24 pm25 pm26 "
            "det1tdc0 det1tdc1 det1tdc2 det1tdc3 det2tdc0 det2tdc1 det2tdc2 det2tdc3 "
            "tc latch\n");
  }

  if (spec_fn) {
    if (!(out_spec = fopen(spec_fn, "w"))) {
      printf("error: could not open output file \"%s\"\n\n", spec_fn);
      exit(1);
    }
    fprintf(out_spec, "# file created by dkl2spec from %s\n", in_fn);
    fprintf(out_spec, "# channel pm11 pm12 pm13 pm14 pm15 pm16 pm21 pm22 pm23 pm24 pm25 pm26 sum1 sum2 sum\n");

    // clear spectra
    memset(pm,   0, sizeof(pm));
    memset(sum1, 0, sizeof(sum1));
    memset(sum2, 0, sizeof(sum2));
    memset(sum,  0, sizeof(sum));
  }

#ifdef HAVE_ROOT
  if (root_fn) {
    // Daniel: first open file, then define trees!
    if ((out_root = new TFile(root_fn, "recreate")) == 0)
    {
      printf("error: could not open ROOT-file \"%s\"\n\n", root_fn);
      exit(1);
    }
    out_root->SetCompressionLevel(compressionLevel);

    // enable automatic splitting of files at <maxTreeSize> bytes
    // TTree::SetMaxTreeSize(maxTreeSize);

    // definiere Tree
    tree = new TTree("Data", "Data");
    // tree->SetAutoSave(autoSave);
    // tree->AutoSave("Overwrite");

    // definiere Branches
    tree->Branch("Cycle", &Cycle, "Cycle/i");
    tree->Branch("StatusCode", &StatusCode, "StatusCode/s");
    tree->Branch("Valid", &Valid, "Valid/O");
    // use signed integers for the latch data : bitwise & operations do not
    // always work properly with unsigned values in Root (changed March/April 07)
    tree->Branch("Latch", &Latch, "Latch/I");
    // added by Bastian 29/03/07 - easy cuts for spin flip states
    tree->Branch("SpinFlipState", &SpinFlipState, "SpinFlipState/O");
    tree->Branch("TimerCounter", &TimerCounter, "TimerCounter/i");
    // added by Torsten 07/03/07 - makes easy time cuts possible
    tree->Branch("DeltaTimerCounter", &DeltaTimerCounter, "DeltaTimerCounter/i");
    //tree->Branch("OldDet10", &OldDet10, "OldDet10/S");

    // Branches related to QDCs and TDCs
    tree->Branch("Det1", &(Pmt[0][0]), "Det1[6]/S");
    tree->Branch("Det2", &(Pmt[1][0]), "Det2[6]/S");
    tree->Branch("DetSum", &(DetSum[0]), "Det1Sum/F:Det2Sum/F");

    // TDC info
    char format[80];
    sprintf(format, "Det1Tdc[%i]/i", TDC_STOPS);
    tree->Branch("Det1Tdc", &(Det1Tdc[0]), format);
    //tree->Branch("Det1Tdc", Det1Tdc, "Det1Tdc/i");
    sprintf(format, "Det2Tdc[%i]/i", TDC_STOPS);
    tree->Branch("Det2Tdc", &(Det2Tdc[0]), format);
    //tree->Branch("Det2Tdc", Det2Tdc, "Det2Tdc/i");
    // added by Bastian 2007/04/12 - easy access to Latch info
    tree->Branch("TdcHits", &TdcHits, "TdcHits/i");

    // Trigger info
    // added by Bastian 2007/04/12 - easy access to Latch info
    tree->Branch("PmtTrigger", &(PmtTrigger[0]), "Det1PmtTrigger/s:Det2PmtTrigger/s");
    tree->Branch("DetTriggered", &(DetTriggered[0]), "Det1Triggered/O:Det2Triggered/O");
    tree->Branch("DetTriggeredFirst", &(DetTriggeredFirst[0]), "Det1TriggeredFirst/O:Det2TriggeredFirst/O");

    // some status report...
    cout << "Creating ROOT file " << root_fn << "..." << endl;
    cout << "  maximum tree size: " << tree->GetMaxTreeSize() << endl;
    cout << "  compresion level : " << compressionLevel << endl;
    // cout << "  autosave every   : " << autoSave << " bytes" << endl;
    cout << "  MaxVirtualSize   : " << tree->GetMaxVirtualSize() << endl;
  }
#endif

  // Cycle time:

#ifdef HAVE_ANSI
  printf("progress:\n  [%51s\n\033[1A[","]");   // Fortschrittsbalken
#else
  printf("progress:\n  [%51s","]");   // Fortschrittsbalken
  for (i=0; i<51; i++) printf("\b");
#endif

  Cycle = header.fcID - 1;
  cyclecounter--;

  // FIXME: this is not meaningful?!
  lastValidCycle = header.cycle_count-1;
  if (header.fileOK == 0) {
    lastValidCycle = header.cycle_count-1;
  }

  ev_count = 0;
  i = 0;
  //2007-02-25, AK:
  // former versions of dackel had a bug, header.ev_count did not
  // contain the correct value, so we do not rely on this
  // for(i=0; i<header.ev_count; i++)
  while(!gzeof(gz_in))  {
    if ((i % p) == 0) {
      printf("*");
      fflush(stdout);
    }
    i++;

    ReadEvent(es, header.version, i);
    Latch = GetLatchContent(es, header.version);

    // detect new cycles
    if ((Latch >> 31) != (prevLatch >> 31)) {
      cyclecounter++;  // index starting from zero
      Cycle++;         // absolute number (offset header.fcID)
      consistency = 0;
      if (header.noc > CYCLE_VALID_COUNTER) {
        // only newer file versions support this!
        // consistency should be zero for valid cycles, non-zero otherwise
        consistency = header.getCounter(CYCLE_VALID_COUNTER, cyclecounter);
      }
      // skip invalid cycles from aborted files (old dackel versions)
      if (cyclecounter > lastValidCycle)
        consistency = 1;  // TODO: set a more meaningfull value here!
      // report skipped cycles
      if (consistency)
        cout << "Warning: cycle " << Cycle << " (index " << cyclecounter << ") has errors, skipping" << endl;
      TimerCounter = 0;
    }

    // skip events iff the cycle is invalid
    if (!consistency) {

      // increase event counter
      ev_count++;

      // timer for the detector stop signals
      GetTdcContent(es, DET1_CHANNEL, Det1Tdc, header.version);
      GetTdcContent(es, DET2_CHANNEL, Det2Tdc, header.version);

      TimerCounter = GetTcContent(es, header.version);
      // 07/03/07 calculate the time difference to last event
      if (DeltaTimerCounter < TimerCounter)
        DeltaTimerCounter = TimerCounter - DeltaTimerCounter;
      else
        DeltaTimerCounter = TimerCounter;

      // extract information from latch bits
      SpinFlipState = (((Latch >> 30) & 0x01) != 0);

      // TODO: move this to new funtion (e.g.: GetTdcHitCount
      // for dackel versions after 20070328 these latch bits contain
      // the number of TDC hits (data format 2)
      unsigned hits = ((Latch >> 27) & 0x07);
      if (header.version == 1) {
        if (hits == 0x04) {
          // bit meaning: more than two and not more than 3
          TdcHits = 3;
        } else if (hits == 0x05) {
          // bit meaning: more than 3
          TdcHits = 4;
        } else {
          // One or two hits, no bits set, so we count ourselves
          TdcHits = 0;
          if (Det1Tdc[0]!=0) TdcHits++;
          if (Det1Tdc[1]!=0) TdcHits++;
          if (Det2Tdc[0]!=0) TdcHits++;
          if (Det2Tdc[1]!=0) TdcHits++;
        }
      } else if (header.version == 2) {
        TdcHits = hits;
      } else {
        // TODO: warning, unknown version!!
      }


      // TODO: move this stuff to new function (e.g. GetPmtTriggerCount())
      // number of PMT triggers
      {
        unsigned count;
        unsigned pattern0 = (Latch & 0x3f);
        unsigned pattern1 = ((Latch>>8) & 0x3f);

        PmtTrigger[0] = 0;
        PmtTrigger[1] = 0;

        for (count=0; count < 6; count++) {
          PmtTrigger[0] += (pattern0 & 0x01);
          PmtTrigger[1] += (pattern1 & 0x01);
          pattern0 /= 2;
          pattern1 /= 2;
        }
      }

      // Latch trigger information
      DetTriggeredFirst[0] = ((Latch & 0x00010000) != 0);
      DetTriggeredFirst[1] = ((Latch & 0x00020000) != 0);
      DetTriggered[0] = ((Latch & 0x00100000) != 0);
      DetTriggered[1] = ((Latch & 0x00200000) != 0);

      // Photomultiplier; ziehe Pedestals gleich ab (sind mit Null initialisiert)
      for (int det=0; det<DET_NUMBER; det++)
        for (int pm=0; pm<PMT_PER_DET; pm++)
          Pmt[det][pm]= GetAdcContent(es, det, pm, header.version) - Pedestal[det][pm];


      if (out) {
        fprintf(out, "%d %d %d %d %d %d %d %d %d %d %d %d %u %u %u %u %u %u %u %u %u %u #0x%08x\n",
              Pmt[0][0], Pmt[0][1], Pmt[0][2], Pmt[0][3], Pmt[0][4], Pmt[0][5],
              Pmt[1][0], Pmt[1][1], Pmt[1][2], Pmt[1][3], Pmt[1][4], Pmt[1][5],
              Det1Tdc[0], Det1Tdc[1], Det1Tdc[2], Det1Tdc[3],
              Det2Tdc[0], Det2Tdc[1], Det2Tdc[2], Det2Tdc[3],
              TimerCounter, Latch, Latch);
      }

      if (out_spec) {
        signed bin;
        int k;

#define inrange(l,x,u) ((x)<(l) ? (l) : ((x) > (u) ? (u) : (x)))
#define binrange(x) inrange(0, (x), (spec_size-1))

        // spectra of all PMs
        for (unsigned det=0; det < 2; det++) {
          for (unsigned pmc=0; pmc < 6; pmc++) {
            bin = binrange(Pmt[det][pmc] / binsize + spec_ofs);
            pm[det][pmc][bin]++;
          }
        }

        // sum spectrum detector 1
        for (bin = 0, k=0; k<6; k++)
          bin += Pmt[0][k];
        bin = binrange(bin / (6 * binsize) + spec_ofs);
        sum1[bin]++;

        // sum spectrum detector 2
        for (bin = 0, k=0; k<6; k++)
          bin += Pmt[1][k];
        bin = binrange(bin / (6 * binsize) + spec_ofs);
        sum2[bin]++;

        // sum of both detectors
        for (bin = 0, k=0; k<6; k++)
          bin += Pmt[0][k] + Pmt[1][k];
        bin = binrange(bin / (12 * binsize) + spec_ofs);
        sum[bin]++;
      }

#ifdef HAVE_ROOT
      if (out_root) {
        // fill Tree
        StatusCode = header.modus;
        Valid = 1;    // TODO: must be updated!!!

        // Detektorsumme (DetSum) und Kalibrierung (DetCal)
        for (UInt_t Det = 0; Det < DET_NUMBER; Det++) {
          DetSum[Det] = 0.;
          for (UInt_t Pm = 0; Pm < PMT_PER_DET; Pm++) DetSum[Det] += Pmt[Det][Pm];
                // ACHTUNG: Detektorsumme wird durch Anzahl der PMTs geteilt!
          DetSum[Det] = DetSum[Det] / PMT_PER_DET;
        }

        // schreibe das event in den Tree
        tree->Fill();
      }
#endif

    } // if (consistency)

    // remember time of last event
    DeltaTimerCounter = TimerCounter;

    // remember last latch value
    prevLatch = Latch;

    // OldDet10 = Pmt[0][0];
  }

  gzclose(gz_in);

  if (out) {
    fclose(out);
  }

  if (out_spec) {
    // write histogramms
    for(i = 0; i <= spec_size; i++) {
      fprintf(out_spec, "%i %u %u %u %u %u %u %u %u %u %u %u %u %u %u %u\n", i - spec_ofs,
              pm[0][0][i], pm[0][1][i], pm[0][2][i], pm[0][3][i], pm[0][4][i], pm[0][5][i],
              pm[1][0][i], pm[1][1][i], pm[1][2][i], pm[1][3][i], pm[1][4][i], pm[1][5][i],
              sum1[i], sum2[i], sum[i]);
    }
    fclose(out_spec);
  }

#ifdef HAVE_ROOT
  if (out_root) {
    cout << endl << "Writing Tree! please wait... " << endl;
    tree->Write();
    // schreibe Tree
    // schliesse output-Datei
    out_root = tree->GetCurrentFile();
    out_root->Write();
    out_root->Close();
  }
#endif

#ifdef HAVE_ROOT
  // we're finished. stop timing
  timer.Stop();
  cout << endl << endl << "real time: " << timer.RealTime() << " seconds, CPU time: " << timer.CpuTime() << " seconds" << endl;
#endif

  //2007-02-25, AK:
  // former versions of dackel had a bug, header.ev_count did not
  // contain the correct value, so we give a warning in case  if(ev_count != header.ev_count) {
  if(ev_count != header.ev_count) {
    printf("\n\n  WARNING: number of events in header (%i) differs from number of events read (%i)!\n\n",
           header.ev_count, ev_count);
  }

  if (cyclecounter != header.cycle_count-1) {
    printf("\n\n  WARNING: number of cycles in header (%i) differs from number of cycles read (%i)!\n\n",
          header.cycle_count, cyclecounter+1 );
  }
  return 0;
}



/*---------------------------------------------------------------------------*/
// Schreibe Counter in ASCII-Datei
// Achtung: Es wird an die Datei angehaengt!
int WriteCounters(DklHeader &header)
{
  stringstream string;

  // oeffne Ausgabe-Datei als stream
  ofstream out(counter_fn, ios_base::out | ios_base::app);
  if (out.bad())
  {
    printf("error: could not open counter output file \"%s\"\n\n",counter_fn);
    exit(1);
  }

  // Get Counter Values and put them on screen
  for (u_int32_t i=0; i<header.cycle_count; i++)
  {
    string.str("");	// string leeren
    string << header.fcID+i << "\t";
    for (u_int32_t j=0; j<header.noc; j++) string << header.getCounter(j,i) << "\t";
    string << endl;
    out << string.str();  // schreibe stringstream in Ausgabestream
  }

  out.close();

  return 0;
}
/*---------------------------------------------------------------------------*/
